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We show that stochastic phase-space methods within the truncated 
Wigner approximation can be used to solve non-equilibrium dynamics of 
bosonic atoms in Id traps. We consider systems both with and without 
an optical lattice and address different approximations in stochastic syn- 
thesization of quantum statistical correlations of the initial atomic field 
operator. We also present a numerically efficient projection method for 
analyzing correlation functions of the simulation results. Physical ex- 
amples demonstrate non-equilibrium quantum dynamics of solitons and 
atom number squeezing in optical lattices in which case we, e.g., numer- 
ically track the soliton coordinates and calculate quantum mechanical 
expectation values and uncertainties for the position of the soliton. 

1.1. Introduction 

Atomic systems with enhanced quantum fluctuations can be prepared in 
tightly-confined cigar-shaped atom traps, where the strong transverse con- 
finement suppresses density fluctuations along the radial direction of the 
trap (see, e.g., Ref. pQ). Quantum effects may be further strengthened 
by reducing the kinetic energy of the atoms by means of applying an op- 
tical lattice potential along the axial direction [21 [3] ■ In this article we 
describe the use of truncated Wigner approximation (TWA) (JHSj m dissi- 
pative non-equilibrium quantum dynamics of bosonic atoms in a Id system 
both with and without an optical lattice. The TWA formalism can include 
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a very large number of degrees of freedom in the stochastic representation 
of the atomic held operator and dissipative dynamics emerges from a micro- 
scopic treatment of the unitary quantum evolution, due to energy dissipa- 
tion within the large phase-space without any additional explicit damping 
terms in the Hamiltonian. Quantum and thermal fluctuations of the atoms 
may be included in the initial state and the resulting quantum statistical 
correlations of the initial state may be accurately synthesized for different 
quantum states in the Wigner representation. Some of the approaches pre- 
sented in this article have been used in the studies of a fragmentation of a 
Bose-Einstein condensate (BEC) by turning- up of an optical lattice 0IH], 
dissipative atom transport [TU] , dynamically unstable lattice dynamics [UJ , 
and dark solitons [TU [JJ5] . 

We address different approximations in synthesization of quantum and 
thermal noise in the initial state. We start with a simple uniform sys- 
tem and phonon excitations within the Bogoliubov approximation. These 
are extended to non-uniform systems and situations where the back-action 
of the excited-state correlations on the ground-state atoms is included in 
self-consistent methods. A particular problem of analyzing non-equilibrium 
quantum dynamics in TWA, related to symmetric operator-ordering of the 
Wigner distributions is addressed by providing a numerically practical pro- 
jection technique. 

1.2. Methodology 

1.2.1. Initial state generation in TWA 

In the approach considered in the present case, Id dynamics is unitary and 
quantum and thermal noise enter only through the stochastic initial state 
that generate the fluctuations of TWA dynamics. The dynamical evolution 
is governed by the Gross-Pitaevskii equation (GPE), but ^w(x,t) should 
be considered as a stochastic phase-space representation of the full field 
operator: 



where the interaction strength g\d — 2huj±a, the s-wave scattering length 
a, the total number of atoms N, and V is the trapping potential. Including 
noise only in the initial conditions is by no means necessary, however, since 
depending on the physical problem, we could include, e.g., atom losses via 
collisions and spontaneous emission that would generate also dynamical 
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noise terms for each time step. Since for the unitary evolution all the 
noise is incorporated in the initial state, it is especially important that the 
quantum mechanical correlation functions for the initial state of the atomic 
field operator are synthesized as accurately as practical for each particular 
physical problem. Here we follow the basic formalism of Refs. [7j [8] which 
has also been reviewed in [TT] and further developed in [T2J [T3] . 

1.2.1.1. Uniform system 

The simplest approach is a weakly interacting bosonic gas in a uniform space 
at T = where we use the Bogoliubov approximation to model quantum 
fluctuations of the atoms. In the Bogoliubov theory we calculate the lin- 
earized fluctuations of the ground state (or a stationary GPE solution) in 
which case the back-action of the excited-state atoms on the ground state 
is ignored [14]. We write the decomposition 

4f{x) = 4> { x )^o+i>'(x), (1.2) 

where the total number of ground state atoms N c — (ajdo). The fluctua- 
tion part tp' for the excited states can be written in terms of quasi-particle 
operators a q and d q as 

$(x) = 4= J2( u « & « eiqx ~ itqt - ^*A^ iqx+Kt ) ■ (i.3) 

The normal mode frequencies e q and the quasi-particle amplitudes u q and 
v q can be solved straightforwardly and the expectation value of the number 
of excited-state atoms in the Bogoliubov theory is 

N' = £(k| 2 + \v q \ 2 )n BE (e q ) + £ K| 2 , (1.4) 

1 1 

with (& q a q ) = nsE^q) = [exp (e q /kBT) — 1] _1 denoting the ideal Bose- 
Einstcin distribution. At T — we have nBE(e g ) = 0. 

In order to construct the initial state for the atoms in the TWA evolu- 
tion we replace the quantum field operators (iff, by the classical fields 
(^W)^w) by using complex stochastic variables (a q ,a*,) in the place of 
the quantum operators (a q , a q/ ) in Eq. (|1.3|) . In the Bogoliubov theory the 
operators (a q ,a q ) form a set of ideal harmonic oscillators and at T = 
(a q ,a*,) are obtained by sampling the corresponding Wigner distribution 
function [TS] 

T^(a 9 ,a;) = -exph2|a 9 | 2 ]. (1.5) 
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In this case each unoccupied excitation mode is uncorrelated with Gaussian- 
distributed noise. The expectation value (a*a q )w — \ specifies the width 
of the distribution and represents vacuum noise, resulting from the symmet- 
ric ordering of the expectation values in the Wigner representation. The 
noise is distributed in space according to the plane waves, with a constant 
density. In the absence of any correlations between the modes, the vacuum 
noise in the uniform space could be replaced by uncorrelated Gaussian noise 
on evenly-spaced numerical grid points. However, the simplest modification 
to the Bogoliubov expansion is to fix the total atom number in each stochas- 
tic realization. This introduces long-wavelength correlations between the 
ground-state mode and the excited-state phonon modes, so that already in 
this simple example there exist non-trivial spatial noise correlations [T2"Ill3| . 
For each stochastic realization the number of excited-state atoms satisfies 

N' = £(K| 2 + \v q \ 2 )(a q a q \) + ]T k| 2 . (1.6) 

9 9 

with the ensemble average N' — J2 q \ v q\ 2 ■ I R Eq. (| 1 . 6|) we have transformed 
the symmetric ordering of the Wigner representation to quantum expecta- 
tion values of normally-ordered operators by subtracting (a*a q )w — 5 from 
each mode. The ground-state atom number is then obtained from the fixed 
total atom number N, so that in each stochastic realization N c = N — N' 
and we set ao = ^ N c + 1/2. The ensemble average of the ground-state 
population N c = N - N' . 

At T ^ we replace Eq. (H) by [H] 

W(a q ,a* q ) = ^tanh(^)exp[-2|cg 2 tanh(£ g )] , (1.7) 

where (; q = e q /2ksT. The Wigner function is Gaussian-distributed with the 
width nBE^q) + 1/2. The formula introduces thermal populations of each 
quasi-particle mode and generates more complex spatial noise correlations. 
After the noise generation, the initial state for stochastic evolution may be 
written as 

¥ w (x, 0) = o (x)a o + -J= ^(u,a,e i!I - v* q a q *e^ x ) . (1.8) 

Here ^w(x, 0) is a stochastic representation of the full field operator for 
the atoms. 



September 7, 2010 0:7 World Scientific Review Volume - 9in x 6in ch'twa'ar 



TWA quantum dynamics 5 

1.2.1.2. Trapped gases 

In a harmonic trap, or in a combined harmonic trap and optical lattice, the 
atoms experience a non-uniform potential that introduces spatially-varying 
initial noise distribution even at T = 0. We write the external potential as 
V(x) — moj 2 x 2 /2 + sEr sm 2 (irx/d), where Er = h 2 ir 2 /2md 2 is the lattice 
photon recoil energy and d is the lattice period. The Bogoliubov equations 
now become spatially dependent and need to be solved numerically [14] . In 
Eq. (jl.3l) we replace u q e lqx / ' \J~L — > Uj(x) and v q e lgx /y/~L — > Vj(x) where the 
index j refers to the mode number. In the lowest order approximation the 
quasi-particle mode functions Uj(x) and Vj(x) are obtained in the Bogoli- 
ubov theory. In several cases of interest where the multi-mode structure 
of the excitations become important, the Bogoliubov approximation is in- 
sufficient due to large contribution of the quadratic fluctuation terms, and 
one consequently needs to use a higher-order theory, such as the gapless 
Hartree-Fock-Bogoliubov (HFB) formalism in which case the ground-state 
and the excited-state populations are solved self-consistently. The coupled 
HFB equations for the ground state and excitations read [El [17] 

(C-U c N c \<t> | 2 ) 0o = O (1.9) 

Cuj - U c N c 4>lvj = ejUj, 

Cvj - U c N c <j>fuj = -ejVj . (1.10) 

where Uj (x) and Vj (x) (j > 0) are restricted to the subspace orthogonal to 
O . Here 

t = -^-^2 + V ^ x ) + 2C/ ciV c |0 o | 2 + 2U e n'(x) - fi (1.11) 

and n is the chemical potential. In Eq. (jl.lOp the Bogoliubov approxima- 
tion is obtained by setting U c — gid, U c = 0. In the gapless HFB theory 
(the Gl version in Ref. [17 ), U c = g\d [l + m' {x) / N c 4>^\ and U c = g\d\ 
here, the depleted density n'(x) — (x)ip' (x)} : and the anomalous pair 
correlation m!(x) = (?/>' \x)ip' (x)) introduce back-action of the excitations 
on the ground-state such that Eqs. (|1.9|) and (|1.10|) must be solved itera- 
tively until the solutions converge. In the non-uniform case the number of 
excited-state atoms is given by 

N' = jdxJ2 Ke( £j ) (K(x)| 2 + Mx)| 2 ) + \ Vj (x)\ 2 ] , (1.12) 

3 

and the total atom number may be fixed in each realization as in the uniform 
case [HH3]. 
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1.2.1.3. Quasi- condensate description 

In tightly-confined Id traps, the phase fluctuations may be enhanced com- 
pared to those obtained using the standard Bogoliubov theory. A more 
accurate description can be calculated using quasi-condensate formalism 
that can be particularly important, e.g., to phase kinks [121 113) . In the 
quasi-condensate description we write the field operator as |18] 

*(x,0) = y/n (x) + Sh(x) exp(i0(x)). (1.13) 

The density 5n(x) and phase 6{x) operators may be written in the 
Bogoliubov-type expansion, requiring Sn/no and I^ZA^I to be much less 
than one, where SI is the spacing on the numerical grid on which we calcu- 
late the operators and A8 is the gradient of the phase operator across one 
gridpoint. Thus (for j > 0) 

9(x) ^E,^-^), (1-14) 
8h(x) — ^n (x) J^j \ Snj(x)&j + Sn*(x)a^ , (1-15) 

where 8j (x) = Uj (x) + Vj (x) and Snj (x) = Uj (x) — Vj (x) are given in terms of 
the solutions to the Bogoliubov equations (see the previous section). This 
results in a stochastic Wigner representation (9w{x), 8nw{x)) of phase and 
density operators. The stochastic initial state for the time evolution then 
reads [12 [13] 

^w{x, 0) = aJ no,w(x) + Snw(x) exp(i9w{x)), (1-16) 
where the ground-state density n 0t w{ x ) — (N c + l/2)|0 o (a;)| 2 . 

1.2.1.4. Relaxation 

We may also consider an ideal, non-interacting BEC as an initial state 
for the TWA simulations, but before the actual time evolution, we can 
continuously turn up the nonlinear interactions between the atoms. If the 
process is slow enough and relaxes to the ground state, we may be able to 
produce the stochastic initial state of the interacting system. Although this 
may simplify the calculations, in practise the technique in a closed system 
does not necessarily converge to the correct interacting state [5]. More 
complex models with open systems, kinetic equations and time-dependent 
noise can help the relaxation process at finite temperatures [19j . 
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1.2.2. Wigner representation and symmetric ordering 

The Wigner distribution returns symmetrically-ordered expectation values 
of any stochastic representations of quantum operators. In particular, the 
expectation values of the full multi-mode Wigner fields in the TWA simula- 
tions of the time-dynamics are symmetrically ordered with respect to every 
mode. In general, this can significantly complicate the analysis of the nu- 
merical results when quantum fluctuations are important [8 . A numerically 
practical transformation of the symmetrically-ordered expectation values to 
the normally-ordered expectation values of physical observables can be done 
using projection techniques (as shown Refs. [7J [5]). In the presence of an 
optical lattice, a natural approach is to project the stochastic field to the 
several lowest mode functions of the individual lattice sites. We define the 
amplitude of the ith mode of the site j as 

a i , j (t)= [ dx[<pf\x,t)]** w (x,t), (1.17) 

J j th well 

where is the stochastic field and ip^ is the ith mode function of the well 
j. The integration is performed over the jth site and the site population 
reads 

(Aj) - X>l.Ai> = E Mi a ^)w - 1/2] (1.18) 

i i 

where (• • • ) denotes the normally ordered expectation value of the quan- 
tum operators, and {■ ■ -)w the symmetrically-ordered expectation values 
obtained from the TWA simulations. Fluctuations are calculated using 
analogous transformations. For the on-site fluctuations of the atom num- 
ber in the jth site we obtain 

(Arij) 2 = (h 2 ) - (hj) 2 

= X/ [( a tj a hj a k,j a k,j)w - ( a i,j a i,j)w(at,j a k,j)w ~ <5ifc/4] . (1.19) 
i,k 

Similarly, the relative atom number fluctuations between the sites p and q 
are obtained from 

[A(h p - h q )] 2 = 22 [{{ai tP a itP - a* q a itq ) (a* k p a kyP - a* k q a k ^)) w 

i,k 

~{ a i,p a i,p — a i,q a h<l)^( a k,p a k,p — ^k,q a k,q)w ~ . (1.20) 
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Alternatively, we could have, for instance, written 

(hj) = I dx{&{x)V(x)) = I dx{^*{x)^{x)) w 

dx (k(*)l 2 - k(*)l 2 ) • (i-2i) 

i 

Calculation of then, however, results in double integrals over the sites 
that can be computationally slow when performed over a large number of 
realizations. 

1.3. Applications 
1.3.1. Dark solitons 

Dark solitons have been actively studied in BECs and in nonlinear optics. 
Although there exist numerous studies of classical solitons, the quantum 
properties of dark solitons are much less known. Numerical TWA sim- 
ulations are suitable for the studies of the creation and non-equilibrium 
quantum dynamics of solitons in Id traps. We consider the experimen- 
tal imprinting method [201 121] , where a soliton is generated by applying 
a 'light-sheet potential', of value V<j, to half of the atom cloud, for time r, 
so that in the corresponding classical case the light sheet imprints a phase 
jump of 4> c = V^t/H at x = 0, preparing a dark soliton. Classically the 
imprinted soliton oscillates in a harmonic trap at the frequency u / y/2 |22j 
with the initial velocity v/c = cos(0 c /2), depending on 4> c and the speed 
of sound c. The soliton is stationary (dark) for cf> c = n, with a zero den- 
sity at the kink. Other phase jumps produce moving (grey) solitons, with 
non- vanishing densities at the phase kink. 

In TWA simulations we generate the initial state using the quasi- 
condensate formalism and vary the ground-state depletion N'/N. At T = 
we keep the nonlinearity Ngid fixed, but adjust the ratio gid/N. This is 
tantamount to varying the effective interaction strength 7j nt = mgid/h 2 n. 
We can also study the effects of thermal depletion by varying T. 

In quantum case, soliton trajectories in TWA fluctuate between different 
realizations. Individual stochastic realizations of l^wl 2 in a harmonic trap 
represent possible experimental observations of single runs. In TWA we 
can ensemble average hundreds of stochastic realizations in order to obtain 
quantum statistical correlations of the soliton dynamics. We numerically 
track the position of the kink at different times in individual realizations 
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and calculate the quantum mechanical expectation values for the soliton 
position (x) and its uncertainty Sx = y/ (x 2 ) — (x) 2 . 




Fig. 1.1. Soliton dynamics in a harmonic trap showing (a-b) the Wigner density 
\^vj(x,t)\ 2 for individual stochastic realizations with the same gidN = WOhujl, <j> c = 2, 
and T = for N = 50, 100, in (a) and (b), respectively; (c) The quantum mechanical 
expectation value for the soliton position (x) (solid lines) and its uncertainty Sx (shaded 
regions) for N = 8000, 440, 50 (curves with decreasing amplitudes) with the same non- 
linearity gidN, (fi c = 2, and T = 0. At N — 8000, Sx is negligible. Quantum fluctuations 
increase Sx and soliton damping, and decrease the speed. 



1.3.2. Atom number squeezing 

In this section we consider an example of a TWA calculation of atom num- 
ber squeezing due to turning up of an optical lattice. Unlike in Refs. [7J [8] 
where a BEC fragmentation in TWA was investigated by a lattice with a 
large number of small sites, we simulate a six-site system with considerable 
multi-mode effects within individual sites. Bose-condensed 87 Rb atoms are 
confined to a cigar-shaped optical dipole trap where an optical lattice is 
applied along the axial direction |23j . The lattice potential is slowly turned 
up from s(0) = 48-Er to s(r) = 96Er. The harmonic trap frequency 
is uj = 2tt x 21Hz, the atom number N ~ 5000, and the lattice spacing 
d ~ 5.7/im. The system was also realized in the recent atom number fluc- 
tuation experiment |24) . 

Due to large individual lattice sites the multi-mode structure of the fluc- 
tuations is important, and the atom number fluctuations are evaluated by 
using the projection technique to several modes in each site, as explained 
in the previous section. The Bogoliubov approximation is not accurate 
due to phonon-phonon interactions and the initial state is calculated using 
the HFB method. The spatially non-uniform distribution of quantum and 
thermal fluctuations is clearly seen in Fig. 11.21 The lowest HFB modes 
dominantly occupy the outer regions of the atom cloud with significantly 
enhanced atom number and phase fluctuations in those sites. Such fluctua- 
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tions could not be represented, e.g., by a uniform stochastic noise sampling. 




Fig. 1.2. The numerical solution of the lowest two HFB modes in a six-site optical 
lattice showing (a) ui(x) (dotted) and vi(x) (solid); (b) U2(x) (dotted) and V^{x) (solid) 
at s = 24 and T ~ 5.5nK; (c) Relative atom number squeezing at different lattice 
height between two central nearest-neighbor sites £ pq = [A(h p — n q )] 2 (n p +n q ) / (4n p n q ) . 
The different data sets correspond to temperatures (from top to bottom) T ~ 5.5nK, 
T ~ 4.5nK, T ~ 4.0nK, and T = 0. 



1.4. Comparisons 

Generating quantum noise in TWA in 2d and 3d may have problems in 
terms of heating the atom cloud during time evolution due to rapid coupling 
dynamics between the modes [6] and divergence of physical observables as 
a function of the number of modes, but Id systems have been more robust. 
Although clearly insufficient, e.g., in a Mott-insulator regime and at very 
low atom numbers, TWA has been successful in describing superfluid dy- 
namics in the presence of considerable quantum fluctuations. For instance, 
TWA simulations |10j were qualitatively able to produce the experimentally 
observed damping rate of center-of-mass oscillations of bosonic atomic cloud 
in a very shallow, strongly confined ID optical lattice, corresponding to the 
dissipative atom transport experiments of Ref. [2] ■ 

The accuracy of the initial state noise generation can be a crucial lim- 
itation, especially in simulations involving very short time dynamics. The 
spatial distribution of phonon excitations in trapped systems can result in 
very rapid noise variation where, e.g., phase fluctuations are dominated 
close to the edges of the atom cloud [5]. For dark soliton dynamics, the 
differences in the soliton trajectories between the cases in which the noise 
was generated within the quasi-condensate description and in the Bogoli- 
ubov theory are notable |12j . Evaluating phonon modes in the linearized 
Bogoliubov approximation may also become inaccurate, compared to self- 
consistent HFB methods. 
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Fig. 1.3. Differences between initial noise generation in TWA. Dark soliton dynamics in 
a Id harmonic trap showing (a-b) the Wigner density \^ w (x, t)\ 2 for individual stochastic 
realizations of TWA for <j> c = 2, T = 0, and N = 50 (parameters as explained in ll.3.H . 
In (a) the initial state is generated within the Bogoliubov approximation and in (b) 
using the quasi-condensate formalism; (c) The quantum mechanical expectation values 
for the soliton position (x) (solid lines) and its uncertainty Sx (shaded regions). The 
lighter curve with larger oscillation amplitude corresponds to the Bogoliubov case and 
the darker one the quasi-condensate case. 
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